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Abstract 

We propose a one-dimensional model of spinor bosons with SU(2) symmetry and a two-body 
finite range Gaussian interaction potential. We show that the model is exactly solvable when the 
width of the interaction potential is much smaller compared to the inter-particle separation. This 
model is then solved via the asymptotic Bethe ansatz technique. The ferromagnetic ground state 
energy and chemical potential are derived analytically. We also investigate the effects of a finite 
range potential on the density profiles through local density approximation. Finite range potentials 
are more likely to lead to quasi Bose-Einstein condensation than zero range potentials. 

PACS numbers: 03.75.Ss, 03.75.Hh, 02.30.IK, 05.30.Fk 
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I. INTRODUCTION 



Integrable one-dimensional (ID) models of interacting bosons and fermions with 5- 
function interaction [l-3] have had a tremendous impact on quantum statistical mechanics. 
In particular, recent breakthrough experiments on trapped ultracold bosons and fermions 
atoms confined to ID have provided a better understanding of quantum statistical effects 
and strongly correlated phenomena in quantum many-body systems. These models contain 
two-body zero range potentials which allows the wavefunctions to be written as a superpo- 
sition of plane waves by means of Bethe's hypothesis |4j. This assumption is true based on 
the fact that every particle can move freely without feeling the presence of others when no 
collision takes place. 

However, Caiogero fl showed that certain mode, with long range potentials can a,o be 
solved exactly, though not using Bethe's hypothesis. He first solved the three-body problem 
with a harmonic potential and a g/r 2 potential, and then generalized it to the A-body prob- 
lem to obtain the exact expression for the ground state energy and a class of excited states. 
Sutherland [6] then derived the exact solutions for the ground state energy, pair correlation 
function, low- lying excitations and thermodynamics of the model with g/r 2 potential for 
both fermions and bosons in the thermodynamic limit by employing the asymptotic Bethe 
ansatz (ABA) which uses Bethe's hypothesis in the asymptotic limit. Since then, many 
models with non-local interaction were solved exactly through the ABA method. Among 
them are the isotropic Heisenberg antiferromagnetic chain |7|], the quantum lattice model 
with nrverse S1 ah square poteMia, g to t - J mode! with long range interaction g the 
nonliner Schrodinger model [10[ and so on. 

The main idea of the ABA is that one restricts oneself to the asymptotic region where the 
particles are considered to be sufficiently far apart, such that their influence on neighboring 
particles is negligible Then one has to show by some unspecified method that the 

system is integrable, i.e., that it has a complete set of independent integrals of motion. For 
example, various authors 12] have shown that for g/r 2 potentials, one can find N integrals 
of motion for the N particle system. Once this is done, one can then conclude that the 
wavefunction is non-diffractive and thus asymptotically given by the BA. Since the exact 
scattering data is known, one can then obtain the exact thermodynamics of the system 



131 ] . It should be pointed out that a common misconception is that the ABA is only a 
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low-density approximation, i.e., N/L — > 0. This is not true and in fact it gives the exact 



thermodynamics for systems with finite density in the thermodynamic limit (see for 
explanations). When using the ABA, the low-density limit N/L — > is only reached when 
the width of the interaction potential between neighboring particles become large. However, 
for the purpose of this investigation, we restrict ourselves to a finite density system where 
the width of the interaction potential between particles is small. A physical example of 
systems with such properties are dilute gases, whose inter-particle interactions are almost 
local. 

In this paper, we investigate the ground state of two-component spinor bosons with 
finite range Gaussian interactions in ID. The interaction potential for this system can be 
expressed in terms of the sum of even powered derivatives of a 5-function. It gives rise to 
certain nonlinear behaviour not observed in systems with spin- independent potentials jl^j ]. 
This kind of velocity- or state-dependent potential leads to more versatility in studying spin 
waves, ferromagnetic behaviour and the relation between superfluidity and magnetism in 



low-dimensional many-body systems, as shown in Ref. 15] for two-component 87 Rb atoms 
on a quantum chip. By using a state-dependent dressed potential, spin degrees of freedom 
in two-component spinor bosons are tunable. This technique for controlling non-equilibrium 
spin motion allows one to study quantum coherence in interacting quantum systems, and to 
experimentally explore predictions of the thermodynamic Bethe Ansatz (TBA) in a system 
of two-component spinor bosons. 

We first introduce the model in Section [Til In Section IHIl we show that the Hamiltonian 
for this model is integrable. In Section IIVI we derive the distribution functions for the 
charge and spin degrees of freedom from the ABA equations. The ground state energy and 
thermodynamics are evaluated in Section [V] in the limits where the interaction strength 
between particles is large and the width of the interaction potential is small. In Section IVIl 
we apply the local density approximation to obtain the density profiles for this model. And 
finally in Section IVHl we conclude with a summary of our results. 

II. THE MODEL 

Let us consider N bosons with SU(2) symmetry confined to a ID wire of length L with 
periodic boundary conditions. Here we denote the internal hyperfine spin states as | t) 
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and | i). The interaction potential between adjacent particles is given by a generic non- 
negative function v(xj — xi) that is even in the inter-particle separation, i.e., v(x) = v(—x) 
and vanishes at large enough distances, i.e., lim x ^ l . 00 v(x) = 0. For such a system, the first 
quantized Hamiltonian is given by 

n = -^J2]^ + 2c £^- - *o - - Nil « 

j=l 3 j<l 

where m is the mass of each boson and c characterizes the interaction strength which is the 
same for all possible collisions, i.e., between two | t) bosons, two | I) bosons, or one | t) 
and one | l) boson. The interactions are repulsive when c > and attractive when c < 0. 
The external magnetic field is represented by H, and the total particle number is given by 
N = N-f + N±. For the rest of this paper we use the dimensionless units of h — 2m = 1 for 
convenience. These units are also used in all figures. 
In the case when 

the model can be exactly solved in the region x 1 <C x 2 <C . . . <C x^ where the width of the 
Gaussian potential a is small relative to the inter-particle separation, i.e., — Xi\ ^> a 
or (N/L)a 1 for every % < N. In this limit, all particles scatter non-diffractively. This 
implies that the asymptotic wavefunction can be written as a sum of iV! terms corresponding 
to the permutations P of the set of asymptotic momenta {ki}. Explicitly, the wavefunction 
can be expressed in Bethe ansatz form as 

= J2MP)eW ^k Pj x^ . (3) 

The argument that supports non-diffractive scattering is as follows. Consider the two- 
body problem N = 2 where the particles are far apart, i.e., x 1 <^ x 2 . Since \x 2 — Xi\ ^> a, 
the particles behave as free particles, therefore the wavefunction is a product of plane waves 
with total momentum and energy given by 

P = k 1 + k 2 , E = kf + k\. (4) 

Through the scattering process, the total momentum and energy have to be conserved. This 
yields a new set of momenta which is either (k[, k 2 ) = (k±, k 2 ) or (k[, k 2 ) = (k 2 , k±). 
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For the iV-body problem, we can think of it as a succession of two particles colliding 
and then scattering to the asymptotic region as free particles, where each two-body collision 
gives rise to a permutation of the momenta. A product of transpositions acting on the 
permutation P leads to another permutation P' . Hence, the scattering is non-diffractive for 
any number of particles. When a — > in the fully polarized case, v(x) — > 5(x) which allows 
us to recover the Lieb-Liniger interacting spinless Bose gas 



III. INTEGRABILITY OF THE HAMILTONIAN 



We know that in the limit a — >• 0, the Gaussian function tends to the 5-function. The 
5-function is not a function in the classical sense, and should be treated as a generalized 



function 



16] instead. Notice that if the potential v(x) is an even function, its Fourier 



transform v(k) = J^° oo v(x)e lkx dx is also an even function, i.e., v(k) = v(—k). This implies 
that the Taylor expansion of v(k) in the neighborhood of k = only consists of even powers 
of k as given by 

oo 

v(k) = Y,hk 2n . (5) 

n=0 

Assuming that the potential meets such restrictions, we can take the inverse Fourier trans- 
form to obtain the potential in position space as 

v (x) = —V / b n k 2n e~ ikx dk 

oo 

= J2 a nS {2n \x), (6) 
n=0 

where a n = (—l) n b n . This result is derived from the fact that the 2n-th derivative of the 
^-function can be expressed as 8 <y2n \x) = ^- /_ (— l) n k 2n e~ lkx dk. 

Let us now consider a Gaussian type potential. The Fourier transform of the Gaussian 
function is still a Gaussian function and is given by 



cxp 



i 2 



_V2^ V 2a 2 
The Taylor expansion of the right-hand side of Eq. (JTJ) at /c = is 



expl-^'l. (7, 



«pI-^) = B-^(t) V - < 8 > 
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From Eqs. (E]) and (jSJ), we deduce that 

It seems a little odd at first glance that an analytic function can be written in the 
form of an infinite sum of generalized functions. We emphasize that this equality does 
not hold at isolated points, i.e., we cannot, for instance, say that the equality holds at 
the point xq. But one can convince oneself that it holds whenever we consider v(x) as 
a continuous linear functional that associates every function ip(x) which vanishes outside 
some bounded region and has continuous derivatives of all orders, a real number (v,ip). 
Mathematically, v(x) is considered a functional in the sense that {v,ip) = j R v(x)ip(x)dx 
where the integration is performed over the real line for this instance. One can also check 
the validity of the expansion v(x) in terms of a linear combination of 8^ 2n \x), denoted as 
v$(x), when ^(x) = J2p A ai ... aN (P\Q) exp(i J2f=i ^p 3 x Qj) ls ^ ne Bethe ansatz wavefunction, 
by comparing the expressions of (v, ip) and (vg, ip)- In Appendix El we verify the claim that 
0,^) = (vs,ip)- 

With this expression for the potential v(x) and after verifying that (v,ip) = (v§,i/)), we 
can re-write the Hamiltonian in Eq. (pQ) as 

j=l i j<l n=0 ^ ' 

Following Gutkin's work 19] . we can show that this Hamiltonian is integrable. The boundary 
condition imposed by Eq. (jit)]) (derived in detail in Appendix [Cj is 

d d\ ( d d\ 



dxj . : dxjj ,x i- x i+i \dx j+1 dxj J " ; " 



„ n\ \ 8 I V dx j+ i dxj 

Here the interaction strength C is now a d x d matrix, where d represents the number of 
internal energy levels. More explicitly, C = cl d where I d is a d x d identity matrix. The 
superscripts + and — on the position of the (j + l)-th particle Xj have the meaning that 
is infinitesimally greater (or smaller) than Xj + \. This boundary condition is a specific 
case of the ones derived in Refs. [l7|, |l8| for velocity dependent 5-function potentials. To 
compute the matching coefficients A(\, fi) and B(X, jj) that are found in Ref. jl9|, we assume 
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that the wavefunctions before collision and after collision are 



M = e^+w+O (12) 



"i+i 



iP\ x . =xf = A(A,/i)e i(A ^ + ^ +l) +5(A,/i)e i( ^ +A ^ +l) . (13) 



Next, we substitute these wavefunctions into Eq. ffTTj) and use Proposition 1 in Ref. |l9| . 
i.e., A(A,/i) + B(\,{i) = 1, which states that there are only two possible plane wave so- 
lutions after collision. These are either where (i) the momenta of scattering particles are 
interchanged, or (ii) the momenta of scattering particles are left unchanged, with the sum 
of their probabilities equal to 1. This yields the solutions for A(X, fi) and B(X, /z), i.e., 

(A-/i)-iCE~o^(-f) n (A-^ 
A(A, M ) = x > ( 14 ) 

A — jl 

B(X,fi) = \ 1 . (15) 

A — jJL 



From Theorem 2(b) in Ref. [19[, the symmetric Bethe ansatz, i.e., Bethe's hypothesis for 
a system of bosons, is satisfied since we have found a pair of commuting matching coefficients 
A(X,fi) and B(X,fi) for any matrix C = cl d . Hence we have shown that this model is BA 
integrable. The N particle symmetric wavefunction can then be expressed as 

ip(x Ql < x Q2 < . . . < x Qn ) = A ai ,„ aN (P\Q) exp li^kpjXQj J . (16) 

P \ j=l J 

This wavefunction is a superposition of plane waves with different amplitudes A ai ... aN {P\Q) 
(not to be confused with the coefficient A(A,/x)) where P and Q are permutations of the 
set of integers {1, 2, . . . , iV}. Each plane wave is characterized by the permutation P of 
wavenumbers {kj}, therefore the sum contains iV! terms. Here a/s represent the spin coor- 
dinates. 

It should be noted that the simple procedure of replacing an analytic function by a linear 
combination of 2n-th order derivatives of the 5-function may lead one to think that any 
Hamiltonian with a pairwise interaction potential which is an even function can be exactly 
solved via the ABA. However, this is not true. The BA integrability conditions met by the 
Gaussian function is actually quite restrictive. First of all, any non-local potential we choose 
has to be well-behaved, smooth and an even function. Secondly, it has to vanish quickly as 
a function of the distance between neighboring particles in order for us to make use of the 



ABA. Thirdly, the Gaussian function is unique in the sense that it satisfies both previous 
conditions, and can still be reduced to a 5-function as its width vanishes to zero. This third 
point enables us to make sure our results reduce to the Lieb-Liniger case in the limit a — > 0, 
which is a necessary condition. These three points eliminate many candidates for a choice of 
pairwise interaction potential. In Appendix [Bj we show that for the case where T = 0, there 
exists a unique solution for the Bethe roots, and that they are good quantum numbers. 



IV. THE GROUND STATE 



The scattering matrix and the ABA equations for this model are derived in Appendix [C] 
and Appendix [D] The ABA equations are given by 

AJ- kj - hi - ic'(kj - ki) n /,-, _ A, 



N 

n 

1=1 



X j- k + \c'(\j - ki) 
A ?: - h 



M 

n 



\i - Xj + ic'(\j - Xj) 
^ Xi — Xj — ic'(Aj — Xj) 



i = l M. 



(18) 



where the effective interaction strength c'{u) = ce~ a u ^ 8 is given in Eq. (1D13I) . Here, M 
denotes the number of spin-down bosons in a system where the vacuum state (initial reference 
state) consists of spin-up bosons. The rapidities for the spin degrees of freedom are given 
by {A,}. 

When T = there are no strings involved in the solution for {Xi}, i.e., all Aj's are real. 
Taking the logarithm of the ABA equations gives 



1=1 



M 



i=i 



dikj — Xi) 
kj Xi 



, (19) 



1 N ( 

9 5^ ^ ( r> 



Xi - k t 



C(Xi - h) 



N 



c'(Xj - h) 
Xi - h 



where 6(x) = 2tan _1 x. Here, quantum numbers Ij are integers (half-odd integers) when 
N — M/2 is odd (even) and J« are integers (half-odd integers) when A/2 — M is odd (even). 
Let us then define the functions h{k) and j(A) to represent "particles" when Lh(k) = 2tiI 
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and when Lj(X) = 2nJ . This yields 

N 



h{k) 



k +TY, 9 (- 



k-k. 



(k - k d ) 



2L^ \c'(k-Xi 



J (A) 



-±5>^1 + 

1=1 " 



c'(k - AQ 



(21) 



c'(A - h) 



1 M / 



A-A, 



c'(A - A,) 



i=i 



c'(A - fe t ) 
A -fa 



n 2 



(22) 



In the thermodynamic limit, 



h(k) = k+ 9 



k-k' 



d(k - k') 



- p(k')dk' 



k-X 
c'(k - X) 



a(X)dX 



- / lnWl + 



d(k - X) 
k-X 



a(X)dX, 



(23) 



X-k 
c'(A - k) 



p{k)dk - / 9 



X-X' 
c'{X - A') 



<r(A')dA' 



- / lnWl + 



d{X - k) 
X-k 



1 2 



p(k)dk, 



(24) 



where p(k) and cr(A) are the distribution functions for charge and spin degrees of freedom, 
respectively There are no "holes" in the ground state, therefore we can safely take p h (k) = 
a h {X) = 0. Define £h(k) = 2ixp(k) and Jrj(A) = 2na(\). Taking the derivatives of Eqs. flU 
and (jUJ) finally leads to expressions for the distribution functions in the form 

p{k) = i- + y K^{k-k')p{k!)dk' -i ^ - A)<r(A)dA + ^ if 2 (Jfe - A)o-(A)dA,(25) 
o-(A) = iy KtiX-fypifydk- J K^X - X')a{X')dX' + J K 2 {X - k)p{k)dk. (26) 

The functions i^i(x) and i^2^) are given by 

^ 2 l 

(27) 



KAx) 



K,(x) 



lc'(x)[l + fx 2 } 

71 [c'(x)] 2 + X 2 

1 c'(x)c'(x)[l + ^x 2 ] = c^x) 
2^~ [c'(x)] 2 + x 2 ~ 2a; 



(28) 
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V. THE THERMODYNAMICS IN THE LIMITS c > 1 AND a < 1 



The model described by the Hamiltonian in Eq. ([I]) does not include any explicit spin- 
dependent forces. Therefore the ground state is ferromagnetic according to a theorem given 
by Eisenberg and Lieb ^0|. When the external magnetic field H > 0, the ground state is 
fully populated by | f ) states which were the reference states that we used to derive the ABA 
equations. When H < 0, all | "f) states will flip into | l) states. The ferromagnetic behavior 
and thermodynamics of the special case a = has been studied in literature 2l|, \22 \. 

When T = and H > 0, our model reduces to the single component case. Here c(A) = 
since the distribution of | \) is zero. Therefore we only have one equation to solve 

p{k) = y^]_ q n wik - k>w + (fc - kr p{k )dk ' (29) 

where ±Q are the "Fermi" points. In FIG. [Hand FIG. [21 we plot p(k) versus k for different 
values of c and a by numerically solving Eq. f[29]) . In both figures, we consider values of a 
and c that are beyond the ABA regime, i.e., values that are outside the limits a < 1 and 
c ^> 1. This is done so that we can more easily visualize how the distribution function p(k) 
varies as both parameters vary. We stress that the curves in FIG. [Hand FIG. [2] become less 
accurate as a tends to larger values or as c tends to smaller values. It is clear from the figures 
that as the interaction width a increases, the distribution of quasimomenta k become more 
centered around the origin. This is because the increase in overlap between single particle 
wavefunctions causes the system to behave more and more like a Bose-Einstein condensate 
where the quasimomenta of particles occupy a smaller region in momentum space. 

Using the relations n = J®q p(k)dk and E/L = J^k 2 p(k)dk, we can approximate p(k) 
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FIG. 1: (Color online) Plots of p(k) versus k for different values of c, with fixed density n = 1. 
The top graph has a value of a = (where one recovers the Lieb-Liniger Bose gas) and the bottom 
graph has a value of a = 0.62. All curves are obtained by numerically solving Eq. (|29p . 
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FIG. 2: (Color online) Numerical plots of p(k) versus k for different values of a, with fixed density 
n = l. The top graph has a value of c = 1.75 and the bottom graph has a value of c = 30. All 
curves are obtained by numerically solving Eq. (|29p . 
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by using Taylor's expansion to get 



p{k) 



2tt 
1 

2^ 



■p(k')dk' 



q it C 2e-^(k~k')y4 + {k _ k ,y 



5a 4 



128ttc J„ Q 



rQ 

/ k l4 p{k')dk' + ... 
J-Q 



1 n 
2n tic 



3a 



5a 4 



1 + "-f k2+ m ki ) + s.Lc 



3Ea 2 



1 + 



5a 2 



5^,4 



k + 



128tt 2 c 



Q 



1 + 



k' 2 p(k')dk' 



+0 



1 



0(a 6 ) 



(30) 



The expression 

I Q Q k' 4 p(k')dk' was evaluated by substituting the dominant terms in p(k') 
into the integral, which gave 



rQ 

J-Q 



k A p{k)dk 



1 



n 



2tx tic 



Q 5 



k [ - — I ) dk = \ 1 -\ 



57T 



2n\ 



c J 



(31) 



To find an expression for the Fermi point Q, we evaluate the integral 

•Q 



n 



p(k)dk 



Q 

7T 



, 2n 3Ea 2 nQ' 2 a 2 Q 5 a 4 

H 1 h 

c 4Lc 



4c 



64nc 



! 2n\ 5EQ 2 a 4 nQ 4 q 4 
+ T / + 32Lc + 64c 



Hence 



Q = nn 



2n 



2n\ ir 2 n 3 a 2 



7r 4 n 5 a 4 
32c 



14n 



4c 

57r 2 n 2 Ea 4 
32Lc 



c / 4Lc 
44n 



3E<\ 2 ( _ 4n _ 3ffa 2 
c 4Lc 



5c 



° (?) + ° ^ ■ (32) 



The ground state energy per unit length of the system is given by 
k 2 p{k)dk 



E 
~L 



Q 



3tt 



2n Q 3 a 2 



1 + — + 



C 47TC 



1 + — ] 



9nQ 2 a 2 Q«a A 7Q 5 a 4 



20c 



167T 2 C 2 647TC 



1 + 



106ra\ 
35c y 



15nQ 4 a 4 " 
448c 



+ 0(^)+0(a 6 ). 



(33) 
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FIG. 3: (Color online) Plot of the ground state energy per unit length E/L versus the interaction 
width a and the interaction strength c for a fixed density n = 1. The surface is generated by 
numerically solving the equation E/L = f®qk 2 p(k)dk. 



Substituting Q into E/L and collecting similar terms yields 
E 1 



-ttV 



] 4 / _ 3\ _ 47i 2 n 2 a 2 f _ 10\ _ 37r 4 n 4 q 4 / _21_ 
"7 \ 7J ^ V "7/ 287 V + 107 



+0(7)+ oK) 



(34) 



where 7 = c/n. With this expression for -E/L, the "Fermi" points can be written explicitly 

as 



Q = Tin 



1-- 1 

7 V 7 



2\ 7r 2 n 2 a 2 



27 



1 



7 



7r 4 n 4 a 4 / 82 \ 

127 v 1 ~ 57; 



1 



(35) 



With the expression for the ground state energy, the chemical potential can be derived 
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a = 0.1279 




FIG. 4: (Color online) Comparison between the analytical results and the numerical results for the 
ground state energy per unit length E/L versus c with a = 0.1279 and n = 1. 

using the relation 

d (E 



dn \ L 



2 2 

7r n 



\ 16 (i _ 15 ^_ 8vr 2 n 2 a 2 / 35 \ _ 27r 4 n 4 a 4 / 189 



37 \ 47/ 57 \ 37/ 77 \ 8O7 



+0 (^) + O (a 6 ) . (36) 

The ground state energy is also calculated numerically for different values of a and c by 
using p(k) in Eq. ( 1291) and the definition E/L = J Q Q k 2 p(k)dk. We thus show a plot of E/L 
versus a and c in FIG. [3j As c tends to infinity, the ground state energy will approach 
7r 2 n 3 /3 as predicted by our analytical results. In FIG. HI we compare our analytical solution 
given in Eq. (I3~4"j) with the numerical solution for the ground state energy per unit length 
E/L when a = 0.1279 and n — 1. It is clear they both agree well when c is large. 
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FIG. 5: (Color online) Axial density profiles from the local density approximation for different 
values of a. Here 7 = 10, total particle number N = 1000, and the density at the center of the 
trap is taken to be n(0) = 1. 

VI. LOCAL DENSITY APPROXIMATION 

In this section, we explore the axial density when the system is confined by an external 
harmonic trapping potential. So far our application of the ABA to solve this model has been 
limited to the case where there is no external confinement. When an external confinement is 
applied, the model is no longer exactly solvable. However, if the external trapping potential 
varies slowly enough, the local density approximation (LDA) [23] can be applied to analyze 
the density profiles in a harmonic trap. 

In the LDA, the chemical potential varies along the axial direction x according to the 
equation 




(37) 
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Using the result in Eq. (|36|) . we then have 



Mo) 



2 2 

moo x 



% 2 n(xy 



1 - 



27T 4 n(x) 4 a; 4 



77 



-fl 
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4 7 


■('* 


189 \ 


8O7/ 



l7r 2 n(x) 2 a 2 ( 35 
57 I 37 



Solving this equation for n{x) gives 



n[x) = n(0)\ 1 



x- 



+ 



7r 4 n(0) 4 a 4 



77 



1 + 



57 

28051 \ 
12007 J 



1 - 
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3t~ 


x 2 " 


)1 







X 



where 



and 



n(0) 



/i(O) 



TV 1-^(1- g)' 



i? 2 



2/i(0) 



To obtain the density profiles, we solve the integral 



n(x)dx = N 



R 



(38) 



(39) 



(40) 



(41) 



(42) 



numerically with total number iV = 1000 particles and particle density n(0) = 1 at the 
center of the trap. 

In FIG. El we show the axial density profiles for different values of a. As the interaction 
width a increases, the particles become more concentrated at the center of the trap in a way 
analogous to a Bose-Einstein condensate. 

VII. CONCLUSION 



In this paper, we studied a system of interacting SU(2) spinor bosons in one-dimension 



with finite range Gaussian potential. Using Gutkin's argument 



19 1, this model is shown 



to be exactly solvable. We applied the asymptotic Bethe ansatz to solve this model when 



the interaction width a is much smaller than the inter-particle separation |xj 



The 



Bethe ansatz equations were derived in Eqs. ( IT71) and (ITBI through the quantum inverse 
scattering method. We went on to derive the particle distribution functions for the charge 
and spin degrees of freedom in Eqs. fl25|) and f )26|) . In the limits c>l, a<l and H > 0, 
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we derived the ground state energy and chemical potential (1361) for the system. The 
spin independent interaction leads to a ferromagnetic ground state. Our analytical results 
were shown to be consistent with the exact numerical results from the asymptotic Bethe 
ansatz equations. Finally, we applied the local density approximation to analyze the density 
profiles of the system in an harmonic trapping potential. From our results, we showed that 
an increase in interaction width a causes the spatial and momentum density profiles of the 
system to more closely resemble that of a Bose-Einstein condensate, in the sense that density 
profiles are more concentrated around the origin. 

This work has been partially supported by the Australian Research Council. 

Appendix A: Proof of (v,ip) = (v$,ip) 



Given the Bethe ansatz wavefunction tp{x) = A (Tl ,,, <TN {P\Q) exp(i J2f=i kp 3 x Q 3 )i ft is 
straightforward to show that 

N 



oo \/2-na 2 



e « 



x 2 Q j2a 2 



^A CTl ... ffJV (P|Q)exp i^k P] x Q] dx Qi 



N 



IV2 



3=1 

1 /' .C 

: exp 



7TQ! 



a 2 k P . 



2a' 



and 



^A (Tl ... (JJV (P|g)exp \ iYl kp i x QiJ exp ( 



/oo 00 , 2\ n ( N \ 

(yj S^(x Qi )^2A ai ... aN (P\Q)ex P 

j+i ) Ln=0 ' ^ ' 

/ N \ 00 ✓ 2 \ ™ 

]>> CT ,..^(p|Q)ex P ij>^ E^r y 

P \ j^i J n=0 ' ^ ' 

/ " \ f a 2 k] 

J2 A at...a N (P\Q)e^p \^Yl kp i X ^J eXP y 2 



(Al) 



Qi 



(A2) 



which verifies the claim that (v,tp) = (vs,ip). 
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Appendix B: Yang-Yang variational principle 



Let us focus on repulsive potentials such that v(x) is positive definite and c > 0. When 
T = 0, Eqs. (USD and flU} reduce to 

N . ( h.-h, \ 

(Bl) 



where Ij is an integer and v(k) is the Fourier transform of v(x). This is the fundamental 
equation for the Bethe roots which can be posed as a variational principle as shown by 
Yang and Yang for spinless bosons 25[]. In order to show that Eq. (IBlj) can be uniquely 
parameterized, we introduce the action 

r N 



with 



$(x) 



2 tan" 1 



x 



CV[X' 



3<l 



dx'. 



Then we need to show that Eq. (IB1I) is given by the minima condition 

dB(ki,...,k N ) 



0. 



To prove this, we further introduce the N x N matrix 



B 



d 2 B 

dkjdki 



"dykj k m ) 



Cp"V 2 {kj km) + {kj km)' 



-2c- 



m - k 



C 2 V 2 (kj - kl) + (kj - h) 2 

which is always positive provided that 



d(k) = v(k) - kv'(k) > 0. 



If that is the case 



lj l Kj 



- h) 



^ c 2 v 2 {k 3 - k{) + {k 3 - hf y 3 



Uj - ui) > 0, 



(B2) 



(B3) 



(B4) 



(B5) 



(B6) 



(B7) 



for arbitrary real {uj}. Hence, the solutions of the fundamental equation exist and can be 
uniquely parameterized by a set of integer or half-integer numbers Ij, as long as = 
v(k) - kv'(k) > 0. 
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We shall exclusively consider such type of potentials. Then, the Bethe roots are real 



numbers from Theorem I on p. 11 of Ref. 24|. Finally if l\ > I m then k\ > k m and if l\ = I m 
then ki = k m as long as tan -1 (k/cv(k)) increases monotonically with k. For the Gaussian 
potential, v(k) = exp(— a 2 k 2 /2), which gives $(k) = v(k) — kv'(k) = v(k)(l + a 2 k 2 ) > for 
all real k. Therefore, there is a unique solution for the BA equations when the Gaussian 
potential is used. 



Appendix C: Derivation of the Scattering Matrix 

We employ the coordinate BA to obtain the scattering matrix between two particles. 
This technique is well known, as used by Yang in solving the spin-1/2 fermion model. 
First consider the region 

R : < x Ql < . . . < x Qj < x Qj+1 < . . . < x Qn < L. (CI) 

Define a wavefunction in R as 

^( x ) = *52 A <n-°N( P \Q) expi(k Pl x Ql + ... + k Pj x Q . + k P]+1 x Q]+1 + . . . + kp N xq N ), (C2) 
p 

where cx^s represent the spin coordinates. This wavefunction is a superposition of plane waves 
with different amplitudes A ai _ aN (P\Q) where P and Q are permutations of the set of integers 
{1, 2, . . . , N}. Each plane wave is characterized by the permutation P of wavenumbers {kj}, 
therefore the sum contains iV! terms. 

Consider a new region R' where particles at position XQ j and XQ j+1 are interchanged, i.e., 

R' : < x Ql < . . . < x Q . +1 < x Q . < . . . < x Qn < L. (C3) 

In this region, the wavefunction is defined as 

^'( x ) = A °i-°n( P \Q') ex P K k Pi x Qx + ■■■ + k Pj x Q]+1 + k P]+1 x Q] + ... + kp N x QN ). (C4) 
p 

From the condition that the wavefunction has to be continuous when xq 3 — > xq j+1 , we have 
the relation 

A ai ... UN (P\Q) + A 

(J1...(Tjv 

(P'\Q) = A (P\Q')+A 

(Tl...(JJV 

(P'\Q'). (C5) 
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where P' and Q' represent the permutations P' = (j j + l)P and Q' = (j j + 1)Q, i.e., 
only the positions of the j-ih and (j + l)-th terms are transposed to get P' from P, and Q' 
from Q. 

The 5- function potential gives rise to a jump in the first derivative of the wavefunction 
at position XQ j = xq j+1 - This jump can be evaluated by considering the Hamiltonian in the 
center of mass frame. In this frame, the new coordinates X and Y are related to the original 
coordinates Xj and Xj+i by the transformation relations 

X=*d^Z±I, Y = x J+1 - Xj , (C6) 

and 

x,=X-^, x j+1 = X + j. (C7) 

Their derivatives are related by 

! <! \ () () 



and 

g d d d_ 1 / g g \ 

Higher order derivatives can be similarly expressed in a straightforward manner. 

The time-independent Schrodinger equation Hip = Eip in these new coordinates is then 
given by 

2ox^- 2 W^ + 2c ^^ Y ^(n + -..U(x,F,x') = ^(x,r,x'). (cio) 

n=0 " \ / J 

where the new set of coordinates X, Y and x' replace the old one x. Also, the dimension 
of x' is less than the dimension of x by two, since we replaced those two coordinates by X 
and Y. Integrating this equation with respect to the Y coordinate from — e to e and then 
taking e — > gives 

dip 



dip 

w 



y= 0+ dY 



y. 1_ /a?y d 2n iP 



n\ V 2 ) dY 2n 

i — u n= g 



(Cll) 

Y=0 



where we have repeatedly used integration by-parts to obtain the right hand side of the 
equation. 
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In the new coordinates, the wavefunctions given in Eqs. (1C2|) and (1C4I) are explicitly 
written as 



j>(X, Y, x') = A ai ... UN {P\Q) expi . . + (k Pj + k Pj+1 )X + ±(k Pj+1 - k Pj )Y + 



(C12) 



and 

ij>'(X, Y, x') = A <n~°N (P\Q') expi (• • • + (k Pj + k Pj+1 )X - ±(k Pj+1 - k Pj )Y + . . )j . 

(C13) 

Substituting the wavefunctions defined in Eqs. (IC12j) and (1C 13|) into Eq. flCllj) separately, 
and then adding both equations together yields the relation 

\{k Pj+1 - k P] ) \A ai ... aN {P\Q) - A ai ... aN {P'\Q)] 
+\{k Pj+1 - k Pj ) \A ai ... aN {P\Q') - A ai ... aN (P'\Q')] 



n=0 



1 / a 



2\ ™ 



(k Pj+1 - k Pj 



2n 



[A ai ... aN (P\Q) - A 

(T1...(TJV 

(P'\Q)) . (C14) 



We introduce the transposition operator Tij which transposes the ith and jth spatial coor- 
dinates of the wavefunction, i.e., 



-4. 



a\...<y N 



... Pi . . . Pj . . . | . . . Qj . . . Qi ■ ■ ■ i 



Qi ■ ■ ■ Qj • • • j 



(C15) 



In matrix form, this operator can be written as Pij]^'"^ = i^,^^,^ Yl^j S arjr7 ' r , 
i.e., Tij = Vij for bosons and = — V%j for fermions where Vij is the permutation 
operator. 

Combining this relation together with Eq. (1C5I) transforms Eq. ( 1C14I) to 

i(k Pj+1 - k Pj ) \A ai ... aN {P\Q) - [T, J+1 ]<;;^ ^...^(P'lQ)" 
i 



n=0 



1 / cr 



ft! 



2 (A;p. +1 -A;p.) 



2» 



A CT ,.. ajv (p|g) + iZ-f4K { ...„> N {P'\Q) .(Ci6) 



Rearranging the terms finally gives us an expression which relates the amplitudes of the 
wavefunction before and after collision, i.e., 



A, 



(T1...CTJV 



(P\Q) 



i{k Pj+1 - k Pj )T jjj+1 + cJ2n=o & ' 



! \ 2 



[^k P3+1 -k Pj )] 2n I 



2n 



(C17) 
^...^(P'lQ). 



-1 (T1...0-JV 
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Here I is the identity operator which is included into the relation so that it can be expressed 
in matrix form. The general expression of the scattering matrix is given by the term inside 
the square bracket as 

Y l0 {u)= h3+ 2 2/s , (C18) 

which relates any two amplitudes before and after collision between particles at the ith and 
jth position whereby the change in momentum is u. The sums in Eq. ( 1C17I) are the Taylor 
expansions of the exponential function given in Eq. (IC18j) . 



For this model to be integrable, the scattering matrix Yij(u) has to obey the Yang-Baxter 
relations. To see whether this is true, we shall consider the transposition of two amplitudes 
through different paths. Without any loss of generality, consider going from t4i23(123|Q) to 
^321 (321 1 Q) along the two different paths 

A 123 (123|Q) = [ri, 2 (A;2-A; 1 )] 213 A 21 3(213|Q) 

= [Yi l2 (^2 - ^i)] 213 ^2,3(A;3 - fci)] 23 % 3 i(231|Q) 

= - fci)] 213 [n, 3 (^ - fci)] 281 ^** - h)] 321 A 321 (321\Q), (C19) 

and 

A 123 (123|Q) = [y 2 , 3 (A:3-A:2)] 132 A 1 32(132|g) 

= [Y 2 Ah - h)} 132 [Y l!2 (k 3 - ^)] 312 ^3i2(312|g) 

= [Y 2 , 3 (k 3 - k 2 )} 132 [Y l!2 (k 3 - k^iY^ih - ^)] 321 A3 21 (321|g). (C20) 

Since the outcome of both paths is the same, they must be equal to each other. In general, 
the scattering matrices satisfy the Yang-Baxter relations 

Y a , b (u)Y c4 (v) = Y Cid (v)Y a!b (u), 
Y a>b (u)Y btC (u + v)Y ajb (v) = Y b7C (v)Y aib (u + v)Y b>c (u), 

Y a , b (u)Y b , a (-u) = 1. (C21) 

Appendix D: Derivation of the Bethe Ansatz Equations 
1. The Quantum Inverse Scattering Method 



We will use the quantum inverse scattering method (QISM) 24| to derive the ABA 



equations for this model. On introducing the operator Rij(u) = PijYij(u) where P^j is the 
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permutation matrix, we have the Yang-Baxter equations in terms of Rij(u), i.e., 

Ra,b( U ) R aA U + v)R b)C (v) = R btC (v)R atC (u + v)R a , b (u). (Dl) 

Notice the difference in subscripts between the above equation and the second equation in 
Eq. (1C21|) . The i?-matrices act on the state space of this N particle system Vn = Yl n =i 
i.e., R a ,b(u) acts non-identically on the tensor subspaces V a and V b and identically on the 
rest of the subspaces. 

Using the Lax representation, we introduce the L-operator which acts on the auxiliary 
space and a quantum state space, i.e., L m {u) = R a ,m(u) where a is the auxiliary space 
and m is the quantum state space. In addition, we also introduce the interwining operator 
R(u) = VR(u) where the permutation operator V has the tensor property on operators 
V(A <g> B)V = B ® A. Hence in Lax representation, the Yang-Baxter relation becomes 

R(u - v)L n (u) ® L n (v) = L n (v) <g) L n (u)R(u - v). (D2) 

The next step is to introduce the monodromy matrix T(u) = Ln{u)Ln-i{u) . . ■ Li(u) 
which is the transition matrix through the entire "lattice". In this form, the Yang-Baxter 
relation can be re-written as 

R(u - v)T(u) <g> T(v) = T(v) ® T(u)R{u - v). (D3) 

Lastly we introduce the transfer matrix t(u) = tr a T(u) where the notation tr a implies that 
the trace is taken in the auxiliary space. As a consequence of Eq. (ID3I) . there exists a family 
of commuting transfer matrices t(u), i.e., [r(u), t(v)] = 0. Following the introduction of the 
operators given above, we can proceed with our derivation of the ABA equations. As stated 
earlier, we are interested in the case where this model has periodic boundary conditions, i.e., 

ip(xi, . . . , Xj = 0, . . . , x N ) = V>(xi, • • • , Xj = L, . . . , x N ). (D4) 

For this condition to hold, the wavefunction defined in Eq. (102[) has to satisfy 

A{Pj, P 1 ,...,P N \Q j ,Q 1 ,...,Q N )= exp(i%L)A(P 1 , ...,P N , P^Qu ...,Q N , Qj). (D5) 

As a result, we obtain 

ex. V (ik 3 L)A E {P\Q) 

= R j+lii (k j+1 - kj) . . . R NJ (k N - k j )R 1J (k 1 - %) . . . R^^k^i - kj)A E (P\Q),(D6) 
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where Ae{P\Q) is the initial amplitude before any transposition. We can abbreviate this 
equation as 

n (k 3 )A E (P\Q) = exp(i%L)A B (P|Q), (D7) 

with the definition 

1Zj(kj) = Rj + i t i(kj+i — kj) . . . RN,j(kN — kj)Ri,j(ki — kj) . . . Rj-ij(kj-\ — kj). (D8) 
If we define the monodromy matrix to be 

T N {u) = L N (k N - u) . . . L 2 (k 2 - u)L x {ki - u), (D9) 
the transfer matrix will have the property 

r(u)\ u=k .^TZ j (k j ). (D10) 
Hence the eigenvalues of Eq. ( ID 71) coincide with the eigenvalues of the equation 



r{u)A E {P\Q) = eMikjL)A E {P\Q) 
at the points u = kj for all 1 < j < N. 



(Dll) 



2. The Algebraic Bethe Ansatz 

The i?-matrix for SU(2) is a 4 x 4 matrix given by 



Rij(u) 



ul — ic'(u)Vij 
u + ic'(u) 



I u—ic'(u) 
u+ic' (u) 




V o 









u 


ic'(u) 


u+\c'{u) 


u+ic'(u) 


ic'(u) 


u 


u+ic'(u) 


u+ic'(u) 









\ 




u— ic'(u) 
u+ic'(u) / 



a{u) b(u) 
c(u) d(u) 



where 

c'(u) = ce- a2u2/8 , 

and the matrix representation of the permutation operator is given by 

( I o o o\ 

10 



(D12) 
(D13) 



10 
\0 1 / 



(D14) 
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Similarly, 



IC[U) 



u + ic'(w) 



/ u— ic'(u) 
u+ic' («) 




V o 





ic'(u) 



u+\c'(u) u+ic'(u) 
u ic'(u) 



u+ic'(u) 





\ 






u+ic'(u) 

g U — ic'(tt) 

u+ic' (u) 



(D15) 



/ 



By choosing the basis for spin-up and spin-down states as 



t) 



I) 







/ \ ! 

we can then act each 2x2 block of the i?-matrix on the spin-up basis vector to get 



u — ic'(-u) / 1 



(D16) 




u + \d(u) \ o / ' 



\C\U) 



u + ic'(-u) I i J 



u 



u + ic'(u) \ Q 

Without any loss of generality, we define the vacuum as 

\ o / 

\/l \ / 2 \ / N 

Hence the action of the monodromy matrix on this state is 



(D17) 
(D18) 
(D19) 
(D20) 





(D21) 



T(u)\Q) = L^h-u) 



A(u) B(u) 



C{u) D(u 





(D22) 



Thus the vacuum \Vt) is an eigenvector of A(u), C(u) and D(u) with eigenvalues Ylf=i a (kj ~~ 



u), and rij=i d(kj — u), respectively Meanwhile, B{u) acts as a creation operator for spin- 
downs. 
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Any arbitrary state $(A) can be created in the form of 

$(A) = B(Ai)S(A 2 ) . . . 5(A M )|fi), (D23) 

where M denotes the number of spin-downs in the system. The action of the monodromy 
matrix on this arbitrary state gives 

/ A(u) B(u)\ 

T(/i)$(A)= w w BA 1 BA 2 ...BA M fi, (D24) 

Since the transfer matrix is the trace of the monodromy matrix over the auxiliary space, we 
only need to consider A(ii)B(\ l )B(X 2 ) . . . B(\ M )\tt) and D{p)B(X 1 )B(X 3 ) . . . B(X M )\Q). 

From the Yang-Baxter equation of the form given in Eq. (ID3j) . we obtain the commutation 
relations 

[A(u), A(v)} = 0, [B(u), B(v)} = 0, (D25) 

[C(u), C(v)} = 0, [D(u), D(v)} = 0, (D26) 

A(u)B(v) = U - V - ic ' {u - v) B(v)A(u) + i ^-^-B(u)A(v), (D27) 
u — V u — V 

/x /x v — u — ic'(v — u) , . . . ic'ff — u) , , , , .„ 
D(w)S(u) = -B(v)D(u) + — — —B(u)D(v), (D28) 

where we took a negative factor in the argument of the i?-matrix because the arguments of 
the i?-matrices in Eq. ( 1Q9|) are negative with respect to u. Therefore 

A( fJl )B{\ 1 )B(\ 2 )...B(\ M )\n) 

= f[ H - A - = = ^ TT ^-f' + ^-^j n) + unwanted terms, (D29) 
5 ' At - A, ^jj l -ki-xc\ii- hi) 



and 



£>(//)B(Ai)B(A 2 )...B(A M )|n> 

= -pr fi - Aj + ic'(/i - Aj ^-r ji-h , n) + unwanted terms> (D30) 

/i - Aj [i - ki - id(n - h) 

The sum of the unwanted terms in Eqs. (1D29j) and (1D30I) vanish when there are no poles in 



the eigenvalue of Eq. fIDllj) 



From Eq. ( 1D11I) . we obtain the ABA equations 

»pomo = - n k i -? + ^ " ti n ^"VT""-' . s = i, . ■ ■ (D3D 

fcj- - *i - 1C ; (% - fc,) A± fcj - A, ; 
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N 



Xi - k + \d{\j - k) 
Xi - h 



M 



Xi — Xj + ic'(Aj — Xj) 



n 



n 



Xi — Xj — ic'(Aj — Xj) 



i = 



1,...,M. 



(D32) 



1=1 



Note that here we cannot make a uniform shift for the set {Xi}, i.e., Xi — > Xi — ic/2 for every 
i, because the effective interaction strength c'(u) depends on the quasimomenta {kj} and 
the rapidities {Xi}. 
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